rm(list = ls())

suppressMessages(library(tidyverse))
suppressMessages(library(magrittr))

c(
list.files("/Pub/Users/wangyk/Project_wangyk/Codelib_YK/some_scr", full.names = T, recursive = T, pattern = "\\.R$"),
list.files("/Pub/Users/cuiye/RCodes/UserCode", recursive = T, full.names = T, pattern = "\\.R$")) %>% 
walk(source)

library(Seurat)
source("/Pub/Users/wangyk/Project_wangyk/Codelib_YK/some_cancers/major_test2.r")
# conflicted::conflict_prefer_all("dplyr", quiet = T)

out_home <- "/Pub/Users/wangyk/project/Poroject/F230731001_OV_sc/"
setwd(out_home)
cs <- unique(unlist(plates))
t_sub <- readRDS("out/SC_part/04.tsub2/t_sub_seuratObj.rds")
m1 <- subset(t_sub,BRCA_mut == "BRCA1_mut")
m2 <- subset(t_sub,BRCA_mut == "BRCA2_mut")
m12 <- subset(t_sub,BRCA_mut != "WT")
m12@meta.data %>% head()


library(SCENIC)


dbFiles <- c(
  "https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-500bp-upstream-7species.mc9nr.feather",
  "https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-7species.mc9nr.feather"
)
# mc9nr: Motif collection version 9: 24k motifs
setwd(file.path(out_home,'data'))

for(featherURL in dbFiles){
  download.file(featherURL, destfile=basename(featherURL)) # saved in current dir
}


# dir.create("cisTarget_databases"); setwd("cisTarget_databases") # if needed

# 找到活性细胞亚群之后，需要提取亚群细胞数据，重新标准化、pca、FindNeighbors、FindClusters，定义Clusters，后续SCENIC分析找每个亚群的调控因子
# 找到如果是上皮细胞，针对本项目，参考原文献中使用inferCNV判断其中的肿瘤与正常细胞群体，然后找到其中的关键转录调控因子

### 筛选数据，构造两个表型文件，其中nGene 计算结果与nFeature_RNA结果一致，所以直接选取 重命名，另外一个记作nUMI,存在int文件夹中
cellInfo <- m12@meta.data %>%
  select(nCount_RNA, nFeature_RNA, CellType) %>%
  rename(nGene = 2, nUMI = 1, CellType = 3) %>%
  select(CellType, everything()) 

mkdir('int')

saveRDS(cellInfo, file=file.path(out_home,"int/cellInfo.Rds"))
cellInfo$CellType %>% unique()

colVars <- list(CellType = cs)
# colVars$EpithelialCell <- colVars$EpithelialCell[intersect(names(colVars$EpithelialCell), cellInfo$EpithelialCell)]
saveRDS(colVars, file=file.path(out_home,"int/colVars.Rds"))



org <- "hgnc" # or hgnc, or dmel
dbDir <- file.path("/Pub/Users/wangyk/project/Poroject/F220531002_KIRC_Cuprotosis_Drug", "data/") 
# RcisTarget databases location
myDatasetTitle <- "T cell" # choose a name for your analysis
data(defaultDbNames)
dbs <- defaultDbNames[[org]]

scenicOptions <- initializeScenic(
    org = org, dbDir = dbDir,
    dbs = dbs, datasetTitle = myDatasetTitle, nCores = 20
)


# 此步骤感觉不是很需要，前面已经弄好了
# Modify if needed
scenicOptions@inputDatasetInfo$cellInfo <- "int/cellInfo.Rds"
scenicOptions@inputDatasetInfo$colVars <- "int/colVars.Rds"
# Databases:
# scenicOptions@settings$dbs <- c("mm9-5kb-mc8nr"="mm9-tss-centered-5kb-10species.mc8nr.feather")
# scenicOptions@settings$db_mcVersion <- "v8"

# Save to use at a later time...
saveRDS(scenicOptions, file="int/scenicOptions.Rds") 
saveRDS(m12, file = 'int/m12_seurat.rds')

m12 <- readRDS('int/m12_seurat.rds')
scenicOptions <- read_rds(file = "int/scenicOptions.Rds")
# 共表达网络
# 使用GRNboost，该工具需要输入表达矩阵，原教程推荐使用counts；
# 在前面设定好之后，这一步是按照规则过滤一些基因，具体阈值可参考原教程
exprMat <- as.matrix(m12@assays$RNA@counts)

# exprMat <- read_rds('/Pub/Users/wangyk/project/Poroject/F220531002_KIRC_Cuprotosis_Drug/int/1.2_corrMat.Rds')
genesKept <- geneFiltering(exprMat,
  scenicOptions = scenicOptions,
  minCountsPerGene = 3 * 0.01 * ncol(exprMat),
  minSamples = ncol(exprMat) * .01
)


# # 在进行网络推理之前，检查是否有任何已知的相关基因被过滤掉（如果任何相关基因缺失，请仔细检查过滤器是否合适）：
# interestingGenes <- as.character(unlist(xgene_lst))
# # any missing?
# interestingGenes[which(!interestingGenes %in% genesKept)]
# # character(0)

# 我们现在可以过滤表达矩阵以仅包含这 3372 个基因。该矩阵现在已准备好进行共表达分析。
exprMat_filtered <- exprMat[genesKept, ]
dim(exprMat_filtered)

# Optional: add log (if it is not logged/normalized already)
exprMat_filtered <- log2(exprMat_filtered+1) 

saveRDS(exprMat_filtered, file="int/exprMat_filtered.Rds")

exprMat_filtered <- read_rds("int/exprMat_filtered.Rds") %>% as.matrix()
scenicOptions <- read_rds("int/scenicOptions.Rds") 
library(SCENIC)

# GENIE3 通常需要几个小时（或几天）才能运行。
set.seed(1110)
Genie3_res <- runGenie3(exprMat_filtered, scenicOptions)


# rm(exprMat)

# GENIE3/GRNBoost 可以检测正负关联。为了区分潜在的激活和抑制，我们将目标分为正相关和负相关的目标（即 TF 和潜在目标之间的 Spearman 相关性）。
# 可以在运行runGenie3的时候运行该步骤，前后都可以
# 计算相关性：
exprMat_filtered <- read_rds("int/exprMat_filtered.Rds") %>% as.matrix()
scenicOptions <- read_rds("int/scenicOptions.Rds") 

runCorrelation(exprMat_filtered, scenicOptions)
corrMat <- read_rds(file = file.path(out_home,'int/1.2_corrMat.Rds'))

corrMat[1:3,1:3]
# 一旦 GENIE3/GRNBoost（和相关性）的结果准备好，就可以运行 SCENIC 的其余步骤。


scenicOptions <- readRDS("int/scenicOptions.Rds")
scenicOptions@settings$verbose <- TRUE
scenicOptions@settings$nCores <- 20
scenicOptions@settings$seed <- 1110

exprMat_log <- readRDS('int/exprMat_filtered.Rds')

# For a very quick run: 
# coexMethod=c("top5perTarget")
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"] # For toy run
# save...

scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions) #** Only for toy run!!
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log)

saveRDS(scenicOptions, file = "int/scenicOptions.Rds") # To save status
scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions,
    exprMat = exprMat_log,
    skipBoxplot = T, skipHeatmaps = T, skipTsne = T
)
saveRDS(scenicOptions, file = "int/scenicOptions.Rds")

# 开始可视化
cellInfo <- data.frame(CellType = Idents(m12))

regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC")
regulonAUC <- regulonAUC[onlyNonDuplicatedExtended(rownames(regulonAUC)), ]

regulonActivity_byCellType <- sapply(
  split(rownames(cellInfo), cellInfo$CellType),
  function(cells) rowMeans(getAUC(regulonAUC)[, cells])
)

as.data.frame(regulonActivity_byCellType) %>% write_tsv("int/regulonActivity_byCellType.tsv")

# 热图
regulonActivity_byCellType_Scaled <- t(scale(t(regulonActivity_byCellType), center = T, scale = T))
library(ComplexHeatmap)
col_fun <- circlize::colorRamp2(c(-2, 0, 2), c("#375e94","white" , "#a20534"))
a <- ComplexHeatmap::Heatmap(regulonActivity_byCellType_Scaled,
  name = "Regulon activity",
  show_column_dend = F,
  row_names_gp = gpar(fontsize = 7),
  col = col_fun,
  column_names_max_height = unit(10, "cm")
)
plotout(p = a, od = "int/", name = "Regulon activity heatmap", w = 6, h = 22)


# 二元表达热图
read.delim("int/3.5_AUCellThresholds_Info.tsv") -> AUCellThresholds_Info

regulonAUC_by_cell <- getAUC(regulonAUC) %>% as.data.frame()
m <- map_df(rownames(regulonAUC_by_cell), ~ {
  #  .x = rownames(regulonAUC_by_cell)[1]

  t_value <- AUCellThresholds_Info %>%
    filter(regulon == .x) %>%
    pull(threshold)
  {
    regulonAUC_by_cell %>% .[.x, ] >= t_value
  } %>%
    as.data.frame() -> d
  
  return(d)
},.progress = T)

rownames(m) <- rownames(regulonAUC_by_cell)
m %<>% mutate(across(everything(),~ as.numeric(.x)))
m %>% rownames_to_column('tf') %>% pivot_longer(-tf,names_to = "cell",values_to = 'activity') -> df
cellInfo %>% rownames_to_column('cell') %>% as_tibble() %>% inner_join(df) -> df

df %>% group_by(tf,CellType) %>% summarise(activity_perc = sum(activity)/n()) -> d
d$activity_perc %>% summary()

d %>% slice_max(activity_perc,n = 4)

p <- ggplot(d) +
  geom_tile(aes(CellType, tf, fill = activity_perc)) +
  theme_bw() +
  scale_fill_gradient(
    low = "white", high = "darkseagreen4",
    labels = scales::percent,limits = c(0,1)
  ) +
  theme(
    axis.text.x.bottom = element_text(vjust = .5, hjust = 1, angle = 90),
    axis.text.y.left = element_text(size = 9), legend.position = "bottom",
    axis.title.x.bottom = element_blank()
  ) +
  guides(fill = guide_colorbar(
    frame.colour = "black", ticks.colour = "black", direction = "horizontal",
    barheight = .5,barwidth = 6, title = "Percentage Activity"
  ))

plotout(p = p, od = "int/", name = "activity_binary", w = 4.75, h = 27.2)


# 寻找特定亚群中的特异性调节因子
# 细胞类型特异性调节因子(基于Suo等人在2018年为小鼠细胞图谱提出的Regulon特异性评分(RSS))。
# 可用于多种细胞类型的大分析，以确定细胞类型的特定规则。
# regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC")
rss <- calcRSS(AUC = getAUC(regulonAUC), cellAnnotation = cellInfo[colnames(regulonAUC), "CellType"])
rssPlot <- plotRSS(rss)
plotout(p = rssPlot$plot, od = "int/", name = "rssPlot", w = 5, h = 12)

plotRSS_oneSet(rss, setName = "C6-CD8+ GZMB (Trm)",n = 5) -> a
plotout(p = a, od = "int/", name = "rssPlot_cd8_trm", w = 4, h = 4)

rssPlot$df %>% filter(cellType == "C6-CD8+ GZMB (Trm)") %>% arrange(desc(RSS)) -> dd

source("/Pub/Users/wangyk/Project_wangyk/Codelib_YK/some_scr/Heatmap_manul.R")


a <- Heatmap_manul(
    data_input =regulonAUC_by_cell[as.character(dd$Topic),], color_used = cs, 
    group_infor = cellInfo,
    Colored_OtherInfor = T,
    saveplot = T, output_dir = "int/", 
    var_name = 'rss_RGN_in_cell_type', width_used = 9, height_used = 4,
    heatmap_name = " ", cluster_name = NULL, show_rownames = T, DoWilcox.test = F,
    rownames_fontsize = 12
)

